arXiv: 1507.06012vl [physics.flu-dyn] 21 Jul 2015 


Implicit Large Eddy Simulation 
of a wingtip vortex 
at Rcc = 1.2 • 10® 

Jean-Eloi W. Lombard* David Moxeyj Julien F. A. HoesslerJ 
Sridar Dhandapani^ Mark J. Taylor^ Spencer J. Sherwinll 


Abstract 

In this article we present recent developments in numerical methods for performing a Large Eddy 
Simulation (LES) of the formation and evolution of a wingtip vortex. The development of these vortices in 
the near wake, in combination with the large Reynolds numbers present in these cases, make these types 
of test cases particularly challenging to investigate numerically. We first give an overview of the Spectral 
Vanishing Viscosity-implicit LES (SVV-iLES) solver that is used to perform the simulations, and highlight 
techniques that have been adopted to solve various numerical issues that arise when studying such cases. 
To demonstrate the method’s viability, we present results from numerical simulations of flow over a NACA 
0012 profile wingtip at Rcc = 1.2 • 10® and compare them against experimental data, which is to date the 
highest Reynolds number achieved for a LES that has been correlated with experiments for this test case. 
Our model correlates favorably with experiment, both for the characteristic jetting in the primary vortex 
and pressure distribution on the wing surface. The proposed method is of general interest for the modeling 



of transitioning vortex dominated flows over complex geometries. 


Nomenclature 



c 

= Wing chord 

T = 

Kolmogorov time scale 

h 

= wingspan 

P = 

pressure 

A 

= aspect ratio of the wing 

P = 

far-held pressure 

Rec 

= Reynolds number based on wing chord 

Coo 

inhow velocity 

x,y,z 

= Carterisan coordinates, 

Poo 

density 

Ax 

= local cell size 

Cp 

Pressure coefficient 

u,v,w 

= Cartesian velocity components 

V = 

kinematic viscosity 


= distance from surface in wall units 

a = 

angle of attack 

t 

= time 

k 

mode 

tc 

= convective time associate to chord 

M 

cut-off mode of the SVV hlter 

At 

= time-step 

P 

P = M — 1 polynomial order of the spectral element. 

lo 

= length scale associate to largest eddies 

Ssvv = 

diffusion from the SVV hlter 

r] 

= Kolmogorov length scale 

Q 

SVV kernel 


1 Introduction 

Understanding the development and growth of wingtip vortices over lifting surfaces is an ongoing research topic 
both in academia and industry. From an academic perspective, fundamental open questions remain, such as 
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the possible re-laminarization of the vortex as it is shed from the wing [T] and the origin of meandering mill], 
the low-frequency movement of the vortex core and the evolution of the vortex structure. Vortices shed from 
lifting surfaces pose challenges to model in many an industrial context such as wind turbines, helicopter blades, 
high-lift configuration of aircraft and high-performance automotive industry EH 13 HIS]. Developing a better 
understanding of the near-wake of the vortex, lying within one chord length of the trailing edge of the lifting 
surface, is therefore essential in understanding the complex flow-structure interactions of interest in these 
problems. The far-field properties of these vortices are also a challenge for the aeronautics industry, where 
their persistence imposes strict limits on distances between landing aircraft m- For these reasons we are 
interested in refining modeling methods for investigating the growth of the vortex in the near-field. 

Conceptually, the simplest approach to ensure that the flow physics are accurately simulated is to perform 
a direct numerical simulation (DNS), in which all necessary scales are resolved at a given Reynolds number. 
For cases at even moderately high Re however, this approach is clearly unfeasible. To demonstrate this, let us 
assume the Kolmogorov hypothesis holds for this flow and as a very rough approximation that the length scale 
Iq associated to the largest eddies is of the same order as the chord length Iq « c. The number of grid points 
needed to resolve the Kolmogorov length-scale relates with the Reynolds number as ry ^ , meaning 

that three-dimensional simulation of a uniformly turbulent flow requires a resolution of R^!^ grid points. For 
aeronautical test cases, where Re is typically 0(10®) or O(IO^), we therefore require 0(10^'^) to 0(10^®) grid 
points to resolve the flow. Even accounting for variations in geometry which may permit varying resolution 
throughout the domain, based on the current rate of advancement of high-performance computing (HPC) 
facilities, resolving fully developed three-dimensional flow at high Reynolds number in a timely manner will 
continue to be well out of reach for the foreseeable future. 



Figure 1: Wingtip vortex developing over a NACA 0012 profile with rounded wing cap in a wind tunnel, 
modeling the experimental setup of Chow et aZ.[T| Both wing surface and streamlines are colored by static 
pressure coefficient Cp = 2 ■ {p — Poo)/PooU^ where = 1 and poo = 1- 

Consequently, there has been ongoing development of modeling methods where small turbulent scales are 
not explicitly computed. Traditional Reynolds Averaged Navier Stokes (RANS) methods, alongside more 
recent advanced such as the Reynolds Stress methods m, have been developed to simulate both the complex 
three-dimensional transitioning boundary layer on the wing and the highly curved flow within the vortex. More 
computationally-intensive methods, such as LES and Lattice Boltzman VLES[T1] have also been developed 
or adapted to investigate such flows, in correlating the simulated results to experimental data. The lattice 
Boltzman method has been used in conjunction with a modified k — e two-equation turbulence model as well 
as turbulence wall shear stress model were used to perform a VLES where the walls were the turbulent flow 
at the wall was modelled. 

The key feature of these studies is in their use of reduced equations or turbulence models, all of which require 
parameters to tune their performance. Since the underlying physical processes that dictate the development 
and evolution of vortices is not well understood, it is therefore difficult a priori to determine appropriate 
settings for these models. The aim of this work is therefore to demonstrate how an implicit LES method, 
in which the number of parameters is comparably very small and is used to provide additional stability, can 
successfully be leveraged to obtain accurate comparisons against experimental data. We appreciate that 
there may be different views of the definition of implicit LES. We have adopted the definition of Sagaut m, 
who explicitly refers to SVV as an implicit LES model and states that “using a numerical viscosity with no 
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explicit modeling are all based implicitly on the hypothesis [...] the action of subgrid scales on the resolved 
scales is equivalent to a strictly dissipative action.” The only influence of the sub-grid scales on the resolved 
scales is therefore dissipative. 

Many existing LES codes are based on finite volume, linear finite element or Cartesian grid methods. We 
will instead investigate the viability of the spectral//ip element method, which lies in the class of high-order 
finite element methods. These methods are widely used in academia, as they offer attractive properties such 
as exponential convergence and low dissipation error for sufficiently smooth solutions M- This is significant, 
as faithfully modeling the wingtip vortex in the far wake is of particular interest in many research communities 
and industries. It is implied that accurate modelling of the far held vortex requires precise modelling of the 
vortex onset in the near held using high precision, low dissipation schemes im mi [Tj nzi ng. However, these 
methods have not been widely applied for the type of industrial, high-i?e cases that we consider in this work. 
Such methods are generally perceived as difficult to implement, as a range of specialised preconditioners, 
mesh generation procedures, parallel communication strategies and stabilisation approaches are needed to 
successfully complete a simulation. 

In this paper we present the SVV-iLES formulation, which utilises the spectral vanishing viscosity approach 
to stabilise numerics m- We show how the issues of implementation and mesh generation can be overcome, as 
well as highlight the benefits that these schemes can have for industrial problems, both in terms of resolution 
power relative to existing studies and computational efficiency. To demonstrate the viability and robustness 
of the scheme, we consider the test case presented by Chow et al. [1], in which the flow over a NACA 0012 
wingtip has been investigated with precise experimental measurements. This case has subsequently become a 
benchmark for vortex dominated flows. To this end we perform an SVV-iLES, at the highest Reynolds number 
considered so far for this case, and correlate our results to the experimental data. As in previous numerical 
studies, in order to reduce the computational cost we have chosen to run the computation at a lower Reynolds 
number. In this work we set Re = 1.2 • 10®, as compared to the experiment which uses Re = 4.6 • 10®. 

In the following section we will briefly report the experimental and numerical methods that have been 
developed and evaluated for investigating the nascent wingtip vortex, as well as the key findings of these 
studies. We then discuss the key features of the SVV-iLES method in Sectionand the numerical methodology 
for our simulations. The evaluation of the SVV-iLES, by close comparison with the extensive data from the 
experiment of Chow et and numerical results of Uzun et aZ.|S], is presented in Section [^before we conclude 
in Section 1^ with a brief overview of the key findings. 


2 Literature review 

We begin with a brief presentation of existing work in the investigation of wingtip vortices; other thorough 
reviews of this held can be found in Rossow |g, Spalart [10] and Green & Acosta [20]. Firstly, experimental 
studies are considered and a summary of the types of flow dynamics that exists for these cases is presented. 
We then discuss results obtained by numerical simulations before emphasizing the possible drawbacks of an 
explicit subgrid-scale model in the context of complex flows, such as the wingtip vortex and highlight the 
originality of the method employed in this study. 

2.0.1 Experimental methods 

For the wingtip vortex the Reynolds number is computed from the chord length c as Rcc = cU^jv, where 
[/oo is the free stream velocity and v is the kinematic viscosity. Giuni EUlll] and Giuni & Benard |23| 
investigated the initial formation and development of the wingtip vortex on a NACA 0012 rectangular wing 
with both square and rounded wingtips at Re^ = 7.4 • 10® for angles of attack a = 0°,4° and 12°. They 
reported that for a given fixed angle of attack, a ‘more’ axisymmetric vortex sheds from the rounded wingtip 
with stronger vorticity within its core compared the square-tip wing. They also reported that for the a = 12° 
case the axial velocity excess is higher for the rounded wingtip. This region is surrounded by a region of 
axial velocity deficit corresponding to rolling up of the vorticity sheet. The weaker axial velocity excess for 
the square wingtip is believed to be a consequence of the more numerous secondary vortices generated by the 
square wingtip. 

Defining the center of the primary vortex as the location of the streamwise helicity peak, they investigated 
the meandering of the vortex in the near held. Two distinct modes of behaviour were observed: for the 
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squared wingtip, the meandering decreased with the distance from the trailing edge whereas for the rounded 
wingtip case, the meandering grew quasi-linearly with distance from the trailing edge. This has also been 
reported by Devenport et al. j2] and Giuni et al. |22| . However Jacquin et al. j3] concluded vortex meandering 
to be insensitive to free stream unsteadiness in the wind tunnel. Fabre & Jacquin |241 and Dieterle et 
al. |15] suggested meandering might also be influenced by the destabilization of a vortex by secondary vortex 
of opposite sense. Jacquin goes further by reporting that the various co-operative interactions between pri¬ 
mary and secondary vortices affected the same range of frequencies as meandering. Zuhal m and Zuhal & 
Gahrib [2S] investigated the correlation between number and strength of secondary vortices and the amplitude 
of meandering. McAlister |29] investigated the wingtip vortex of NAGA 0015 with both rounded and square 
wing cap reporting the maximum azimuthal velocity of the vortex to be independent of Reynolds number but 
dependent on the angle of attack. 

Finally, Ghow et al. [I] performed an extensive experimental and numerical study of the wingtip vortex 
over a rounded NAGA 0012 profile in view of defining a benchmark for numerical models. They reported 
merging of the secondary and tertiary vortices into the primary within one chord length of the trailing edge. 
From this distance onwards, the vortex had an axisymmetric structure with a jet-like axial velocity profile 
(i.e. the axial velocity within the vortex core was greater than the freestream velocity). The peak jetting 
velocity was measured to be 1.77t/oo at the trailing edge and by three-quarters of a chord length downstream 
of the trailing edge it was measured to be 1.7Uao. Pressure taps were placed on the wing-surface to assess the 
state of the boundary layer, both under the primary vortex where a strong suction region is present and at 
different spanwise locations, allowing for insight into the three dimensional boundary layer. They also reported 
skin-friction lines to gain understanding into the attachment and detachment lines of the strongest vortices in 
the near wake, as well as low amplitudes of meandering because of both a relatively large angle of attack of 
a = 10° and measurements at relatively short distances from the trailing edge. Devenport et al. |2] reported 
meandering amplitude growing linearly with distance from trailing edge. 

2.0.2 Numerical methods 

Presently, RANS based methods with linear eddy viscosity models, such as fc — w SST and k — e, remain 
commonplace for industrial flow simulation m- For vortex dominated flows, and more generally flows with 
strong curvature, these models may struggle due to the largely unsteady dynamics of the flow. Large discrep¬ 
ancy with experimental data, exceeding 100% error on Cp distribution in the vortex core, as well as drastic 
under-prediction of the jetting within the core have been reported by Ghurchfield et al. [3D]. These results 
have motivated the development, over the past 30 years, of more complex closure models tailored for highly 
curved flows such as the wingtip vortex. We present a brief overview of this development. 

Dacles-Mariani et al. |D] ran a 5*^*^ order compact (i.e. 7 point stencil instead of 11) biased upwind scheme 
for the advection term and second order scheme for the viscous term. Here they underline the necessity 
for low numerical dissipation. They ran a modified version of the one-equation Baldwin-Barth turbulence 
model where they modified the production term to avoid overproduction of eddy viscosity in the vortex core. 
They successfully captured a secondary structure and computed the axial velocity profiles of the core to 
within 3% of experiment but under-predicted the core pressure by more than 25%. It should also be stressed 
that Dacles-Mariani et al. ED used experimental data to setup both inflow and outflow boundary conditions 
for their RANS simulation. Linear eddy-viscosity being too dissipative. Graft et al. [32] developed a non¬ 
linear eddy viscosity model (EVM). This more advanced model still suffered from a more severe decay of the 
turbulent stresses than measured experimentally. Duraisamy and laccarino |33| modified the eddy viscosity 
coefficient of the — f turbulence model and compared their results favorably for axial surplus compared 
to the baseline Spallart-Allmaras and Menter’s k — oj SST models. The wingtip vortex exhibits a peculiarity 
in the turbulence structures where the stress and strain are out of phase which renders questionable the use 
isotropic eddy-viscosity based prediction methods such as (k-w, etc). Ghurchfield et al. |3411351 [T51 111] modified 
the Spalart-Allmaras model to account for streamline curvature and successfully modeled the lag between the 
mean strain rate and respective Reynolds stress. This method produced the best correlation with experiment, 
for a RANS based method, but proved costly (relative to other simpler RANS methods) and so their use may 
be restricted to flows dominated by vortices. They also showed that without the accurate modeling of the 
three dimensional boundary layer the developing vortex remained challenging to compute accurately even for 
the advanced RANS models correcting for the high degree of curvature in the flow. 

In an attempt to develop increasingly robust models, more advanced numerical methods such as Large 
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Eddy Simulation (LES) and Very Large Eddy Simulations (VLES) have been proposed to investigate un¬ 
steady features of the flow as well as aero-acoustic properties. Fleigh et al. [SS] developed a compressible LES 
to investigate far-fleld broadband noise generated by the nascent vortex on a rotating wingtip at Reynolds 
Rcc — 1 • 10® and reported computed power and thrust coeffldents for the windmill blade within 3% of ex¬ 
perimental data. Ghias [S] reported a compressible LES of NACA2415 with square tip run at Rbc = 10® 
where they employed a dynamic sub-grid scale model but did not present the correlation of their data with 
experiment. Uzun et al. ElE] numerically investigated Chow’s experiment with a compact finite differencing 
LES with implicit spatial Altering at Rcc = 5 • 10®. Their results were included in the following comparison 
of our simulations with experimental data. Jiang et al. m reported results for a LES simulation of Chow’s 
experiment at the experimental Reynolds number of Rcc = 4.6 • 10® but did not compare the results with 
experiment. The lattice Boltzman method was used in conjunction with a modified k — e two equation turbu¬ 
lence model and wall shear stress model to perform a VLES. Despite good correlation with the experimental 
results of Chow et al. [T] for the suction of the vortex on the surface of the wing the method over-predicted the 
jetting phenomena within the vortex by 23% at streamwise location x/c = —0.114 and 12% at x/c = 0.456. 
So far these more advanced modeling methods have generally been used for simulations at lower Reynolds 
numbers I371E] or have not not been compared against experimental data in the way they will be in the 
present study 

2.0.3 Motivation and contributions of present study 

Modeling unknown physics by leveraging a sub-grid scale model outside of its operational window can be seen 
as a substantial drawback for explicit sub-grid scale models, which tend to rely on a wide range of parameters 
to dictate their behaviour. The two-equation eddy-viscosity k — uj SST turbulence model introduced by 
Menter [3S] in 1994, for example, relies on 9 modeling constants. Menter underlines, in this paper, the strong 
sensitivity of the resulting computed flow to variations of 5-10% of these constants and further stressing “None 
of the available theoretical tools (dimensional analysis, asymptotic expansion theory, use of direct numerical 
simulations (DNS) data, renormalization group (RNG) theory, rapid distortion theory, etc.) can provide 
constants to that degree of accuracy.” Leveraging these many parameter sub-grid scale models, such as k — uj 
SST and more recently Reynolds Stress Relaxation models m, in complex flow cases therefore requires an a 
priori knowledge of the flow physics. This is often infeasible, particularly when complex geometries are present 
and the length scales of both flow and geometry vary substantially. For example, in engineering applications 
such as flow over a Formula 1 car, different sets of parameters may be required to accurately model the various 
flow dynamics that are induced by the variations in geometry across the body of the the car. This includes 
the vortex dominated flow of the wing-tip regions of a front wing, regions of the front wing where the flow 
is mostly two-dimensional, and the rotating wheel that impinges on the moving road. It may therefore be 
impossible to obtain values for these parameters that capture the desired flow features across the entire car. 

The aim of this paper is to show that regularized high order spectral//ip element methods, without an 
explicit sub-grid model, can be applied to produce results that compare favorably against experimental data, 
by considering flow over a three-dimensional geometry of practical interest. The SVV-iLES approach, which 
we will present in the following section, requires the choice of two regularization parameters: one to dictate the 
level of artificial viscosity, and another for a cut-off wavenumber. These are chosen through experimentation, 
such that the computation does not diverge but do not require assumptions regarding the physics of the flow. 
We discuss this methodology in the following section before showing results of the NACA 0012 wingtip vortex 
case. 


3 Computational Methodology 

In this section, we give a brief summary of the computational methodology used for the NACA 0012 wingtip 
simulations. We begin by outlining the Nektar++ spectral//ip element framework in which the solver is 
implemented |39| . We then outline the types of regularization that are necessary to perform the computations 
and prevent the simulation from diverging, along with the mesh generation procedures that are used to generate 
a curvilinear boundary layer mesh for the geometry. Finally, we discuss initial and boundary conditions as 
well as resolution requirements in space and time for the simulations. 
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3.1 NektarH h: a high-order spectral//ip element framework 

High-order finite element methods often suffer from the stigma of difficulty of implementation, which in turn 
means that despite their attractive numerical properties and the ability to resolve difficult cases such as the one 
presented here, they are frequently under-used. NektarH—h is a framework designed to address this problem 
by providing a modern development environment for these methods. It is highly parallel, providing a range 
of efficient preconditioners and has support for a variety of solvers including the incompressible Navier-Stokes 
equations. For a summary of functionality one can see for example Cantwell et al. [39], or for more details on 
the method itself, the reference book by Karniadakis & Sherwin |40| . 

In the following sections, we outline the modifications made to Nektar++ in order to adapt the existing 
DNS solver, making it suitable for an iLES approach through regularization. Furthermore we outline the 
other challenges that need to be addressed in order to more generally make spectral//ip element methods 
viable for these problems. 

3.2 From DNS to SVV-iLES: filtering with Spectral Vanishing Viscosity 

Running cases at high-Reynolds number on an under-resolved mesh requires close inspection of the source 
of errors. The highly non-linear nature of the underlying equations leads to a complex interactions of these 
errors, which when left uncorrected leads to a diverging solution. Broadly, we have found that the two most 
important aspects are: 

1. consistent integration of non-linear terms; 

2. artificial dissipation to prevent divergence of the flow in the presence of under-resolution. 

We now explain each point in more detail. The incompressible Navier-Stokes solver used in Nektar++ 
directly integrates the underlying equations through the use of an operator splitting scheme in combination 
with a consistent boundary condition for the pressure Poisson equation m- In this scheme, nonlinear terms 
are computed explicitly at each quadrature point, which depending on the element type uses a form of Gauss 
quadrature. These nonlinear terms are then multiplied by the elemental basis functions and integrated in 
order to compute the inner product , as is required in the continuous Galerkin formulation. However, 
since Gauss quadrature will only produce exact values for integrals of polynomials of degree 0{2P) at a 
simulation polynomial order of P, aliasing errors are introduced due to the quadratic nonlinearity present 
in the Navier-Stokes equations. When simulations are adequately resolved, this aliasing error usually does 
not affect the robustness of the simulation. However, when under-resolution is used for implicit LES, this 
aliasing effect leads to a significant buildup of error as the simulation progresses in time, and usually causes 
the simulation to abruptly diverge. 

We note that the nonlinear terms of the Navier-Stokes equations are consistently integrated if the elements 
are straight-sided. In this case, the Jacobian of the mapping which defines the coordinates of the element is 
affine with constant determinant. However, where the element is curved, this mapping is an isoparametric 
polynomial expansion, which when incorporated into integrands, leads to an additional source of aliasing error. 
In this work, we do not take this source of error into account. A more detailed description of these aliasing 
errors, as well as means of suppressing them, can be found in |35|. Whilst it is true we could likely suppress 
aliasing errors using SVV, which we describe below, it would require stronger diffusion together with a lower 
cut-off mode, leading to reduced accuracy of the overall solution. Additionally, since the SVV operator is 
anisotropic, whereas dealiasing is isotropic, the two do not completely overlap. Therefore a reduced amount 
of regularisation can be achieved using dealiasing. 

In regards to the second point, we first note that the energy spectrum of the flow consists in a resolved 
range of wave numbers, the large eddies, and an un-resolved range, the turbulent or dissipative scales, for 
higher wave numbers. Because of this cut-off the higher wavenumber dissipative scales are not resolved for 
LES. The energy build-up at high-wavenumber and the coupling through the nonlinear term down to lower 
wave numbers may lead to unstable computations at worst and erroneous energy spectra at best. 

Spectral Vanishing Viscosity (SVV), first introduced by Tadmor [19] for spectral Fourier methods, aims 
to damp the high-wavenumber oscillations without impeding the physics of the flow at lower wavenumbers, 
thereby stabilising the simulation and preserving the accuracy of the solution. In this approach one adds an 
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additional reaction term of the form 


d 

esvv^ 

ox 


Q * 


du 

dx 


where esvv is a constant, * denotes the convolution operator and Q is a kernel dictating which modes receive 
damping. This approach was extended to the Navier-Stokes equations by Kirby & Sherwin |43| . Karamanos 
and Karniadakis m and has been extensively used by Pasquetti et al. jJS], Severac and Serre [15], Xu et 
al. m as well as Lamballais et al. |48| . 

We also stress that the oscillations stabilized by the SVV method are sub-element oscillations. The intrinsic 
nature of the spectral//ip element methods results in degrees of freedom within each element. When the 
solution field is under-resolved, it is no longer guaranteed to be smooth. This therefore leads to the development 
of spurious high-frequency oscillations in the polynomial representation of the solution inside each element, 
arising from Gibbs phenomena occurring between two connected elements. The aim of SVV stabilisation is to 
prevent these sub-element oscillations, since they may lead to oscillations occurring in neighbouring elements 
and ultimately to divergence of the computed solution. 

In the context of the simulations presented here, we introduce artificial damping as an implicit sub-grid 
scale model through the Spectral Vanishing Viscosity (SVV). We appreciate that there may be different views 
of the definition of implicit LES. We have adopted the definition of Sagaut m, who explicitly refers to SVV 
as an implicit LES model and states that “using a numerical viscosity with no explicit modeling are all based 
implicitly on the hypothesis [...] the action of subgrid scales on the resolved scales is equivalent to a strictly 
dissipative action. ” The only influence of the sub-grid scales on the resolved scales is therefore dissipative. 

The key point in SVV filtering is that, due to the shape of the kernel Q{k), 


Q{k) = 


cxuf 

t^xp (M-k)2 J > 

0 , 


k > M, 
k< M, 


artificial viscosity for any mode number k is only applied above a cut-off mode M. For the higher modes, the 
total viscosity can thus be expressed as l/i?e -I- esvv- 

The SVV operator was incorporated into the velocity correction scheme of the incompressible Navier-Stokes 
equations inside Nektar++ by following the approach presented by Kirby & Sherwin jJS], where the elemental 
Laplacian operator is convolved with the kernel Q. The computational cost is therefore negligible since the 
only cost involved is during setup. However, we note that the addition of SVV can lead to higher iteration 
counts in the conjugate gradient method used to calculate the intermediate pressure field and perform the 
velocity correction. This effect can be mitigated through the use of appropriate preconditioning strategies |49| . 

In our SVV-iLES method the parameters do not adapt automatically to the flow (i.e. without the input 
of an experienced user), although methods have been proposed to overcome this constraint by implementing 
an adaptive SVV diffusion muni- These methods however still require the same number of parameters to be 
calibrated on a case-by-case basis and increase the computational runtime cost, as the matrix systems which 
represent the diffusion operator need to be rebuilt whenever either parameter is changed. We have therefore 
not considered such approaches here. Also note that we do not use spatially-variable SVV diffusion coefficient 
or wall functions to limit the amount of damping near solid walls. 

As reference, the values of the SVV parameters used by different groups, including the simulations per¬ 
formed here, are reported in Table Since we do not have a detailed a priori knowledge of the flow features 
and, consequently, of the sub-grid flow physics, we have chosen our SVV parameters arbitrarily so that they 
ensure a non-diverging solution. For other simpler cases in the literature, more effort has been made to more 
rigorously quantify the use of specific SVV parameters |48| . For the sake of comparison however, we report the 
parameters used in previous studies and stress that our choice aligns reasonably closely with that of previous 
studies. We recognise greater exploration is necessary for more complex cases. Research along these lines is 
presently being conducted by Moura et al. m- 


3.3 Geometry and mesh generation 

The rectangular wing investigated numerically has a NACA 0012 profile, with a rounded wing cap (conse¬ 
quently a longer semi-span where the wing is thickest) and a blunt trailing edge. The semi-span, without 
the cap is, b = 0.91[m] and the chord is c = 1.22 [to] which correspond to an aspect ratio of A = 0.75. The 
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M 

esvv 

TadmorfTI?] 


l/iV2 

Pasquetti|3Sl HU HU 

{Vn, In, ifv, IN} 

{1/A^,1/47V,4/1V} 

Xu[i7| 

{N-2,lN} 

1/N 

Karamanos|H] 


1/N 

KirbyHO] 

bVN 

5/8 

Present study 

0.5N 

0.1 


Table 1: Different values for the SVV parameters (diffusion esvv cut-off mode M) where = P — 1 for 
a discretization of polynomials of order P. 



Figure 2: Computational domain representing the test section in the wind tunnel used by Chow et al. [T] for 
their experiment. 












boundary layer tripping mechanism [T] used in the experiment is not reproduced in the mesh. We represent the 
test section of the low speed wind tunnel located in the Fluid Mechanics Laboratory of NASA Ames Research 
Center used by Chow et al. [1] for the experimental work, by a 0.66c x c x 10c cuboid domain as shown in 
figureWe do not model the wind tunnel sections upstream and downstream of the test section. A Cartesian 
coordinate system {x, y, z) is used to locate point within the computational domain with its origin the wingtip 
trailing edge, where x is aligned with the streamwise direction and y, z are the two transverse ordinates. 

Regarding the interference between the primary vortex and the wind tunnel walls, Chow et al. decided 
to use as large a model as possible while still avoiding severe viscous interaction between the wall and the 
primary vortex through significant growth and/or separation of the boundary layers on the wind tunnel walls. 
The authors of the experiment warn against large inviscid effects (mirror effects) due to the close proximity 
of the wall that most likely influence both the primary and secondary vortices. It should also be noted that 
the significant blockage created by the large wing in the relatively small wind tunnel section accelerates the 
flow around the wing. Hence the absolute angle of attack, perceived by the wing, is around 2° higher than 
the angle of attack prescribed by the geometry. 

Meshing methodology When considering complex geometries, even generating a linear, straight-sided 
mesh poses a significant challenge. We have therefore turned to commercial mesh generators, which provide a 
robust approach to generating linear straight-sided meshes, but generally lack support for high-order elements. 
We note that for high-order simulations, elements which lie on the boundary must be curved so that they align 
with the underlying geometry. Using straight-sided high-order elements can significantly alter the physics that 
form near boundaries and thus downstream of the boundary. 

The commercial software we use first imports the CAD geometry, in this case the wing, and makes a 
fine tessellation of the surface. This tessellated surface is then used to produce the surface mesh. As we 
have the meshed surface and the original ICES (or CAD) geometry, but not the intermediary tessellation, we 
are presently unable to add the necessary curvature to the linear mesh by interrogating the CAD geometry 
directly. To smooth the mesh we adopt an alternative, patch-based technique known as spherigons |54| which 
rely on surface normals that are obtained from a fine triangulation generated by the commercial software. 

The robust high-order three step mesh generation procedure used is summarized in Fig. [^and the interested 
reader should refer to Moxey et al. for further details. Our meshing procedure generates a coarse single¬ 
element boundary layer of prisms which are curved using spherigons at the wing surface, with straight-sided 
tetrahedra filling the rest of the volume. Each prism is split using an isoparametric method in the wall-normal 
direction, which allows us to achieve the desired wall normal resolution whilst preventing the self-intersection 
of the boundary layer elements. The tetrahedra grow from the prism layer in a controlled manner in regions 
of interest, such as the vortex path, where we impose a regular element size to avoid introducing additional 
mesh induced error. 

Comparative degrees of freedom The implicit diffusion from the SVV operator preserves the convergence 
properties of the underlying scheme. It does not degrade the exponential rate of convergence of the accuracy 
of a solution achievable for a sufficiently smooth field m- We have focused on using our anisotropic prism 
refinement technique to capture the wall-normal near-wall resolution. In this region we have measured that 
across the wing surface, the placement of the closest grid point satisfies y~^ < 1. With this level of resolution, 
we can accurately capture the sub-viscous layer of the flow field and resolve the high shear of the boundary layer 
profile. However the resolution of subsequent elements is not enough to, for example, capture the extremely 
fine scale of the turbulent boundary layer characteristics. Additionally, for computational reasons, we clearly 
cannot resolve with a similar level of accuracy in the wall-tangential direction. We do appreciate that this 
may very well play an important role in capturing the flow in the region of the vortex roll-up. We discuss this 
further in the presentation of our results. 

The present mesh is composed of 243,000 elements of which 24,500 are prisms around the wing surface 
and 218,500 are tetrahedra growing from the three prism layers to the wind tunnel walls. The prism layer 
represents roughly 20% of the total number of degrees of freedom. Running this computation with 5*^ order 
polynomials (i.e. 6*^ order accuracy in space) amounts to roughly 16.7 million degrees of freedom; around 
40% fewer degrees of freedom than used by Uzun et al. for their LES. Table compares degrees of freedom 
for our SVV-iLES method with other studies of this case. 
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Figure 3: Robust three-step meshing procedure. 


Method 

DOF (TO®) 

RANS (modified Baldwin-Barthl |55] 

2.5 

RANS (Lag RST)[TI] 

13.8 

LES IS] 

26.2 

iLES |37] 

26 

Present study 

16.7 


Table 2: Comparison of mesh used by different methods, converted to degrees of freedom (DOF) 
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Figure 4: Overview of the coarse surface mesh and two planes, spanwise cut at the root of the wing and 
streamwise cut around mid chord. For sake of clarity quadrature/collocation points within each surface 
element are not shown here. 

3.4 Initial and boundary conditions 

The simulation is impulsively started from u/Uco = 1 throughout the domain, except at the no-slip boundaries 
where u/Uoo = 0, at a low Rcc « 10. We then gradually increase the Reynolds number by a factor of 10, 
each time waiting for two convective time units defined as tc = c/Uao, until the simulation reaches the 
desired Reynolds number. To improve computational efficiency the polynomial order P of the discretization is 
increased with Reynolds number, with P = 2 for Rbc = 10 and P = 6 for Rec = 1.2 • 10®. At the outflow, we 
impose the boundary condition developed by Dong et al. m that balances the kinetic energy influx through 
the outflow boundary condition to prevent instability. The computational setup differs from the experimental 
setup in three ways. We discuss these and their possible influence on the results in the following paragraph. 

Firstly, the boundary layers developing on the wind tunnel walls are neglected, by using a free slip condition. 
The primary vortex is located sufficiently far away so that the viscous interaction between the primary vortex 
and the wind tunnel walls is much weaker than the interaction between the primary vortex and the wing 
surface via the secondary structures. 

As a first approximation, the wall acts, inviscidly, as a symmetry condition. Since the computed location 
of the vortex is similar to the experimental data, the inviscid interaction between the vortex and the wind- 
tunnel walls is assumed to be of similar intensity. Secondly, as with other LES studies, we do not model 
turbulence at inflow. As we shall see in the results section, despite the lack of a turbulent inflow we still see 
good comparisons against experimental data. Finally the boundary layer on the wing is not tripped. The 
tripping of the boundary layer near the leading edge has been reported to increase the diameter of the vortex 
(measured by the peak to peak distance of the vertical velocity) by 30% [1^. McAlister also report that 
adding a boundary layer trip changes the streamwise component of the velocity from a small excess to a large 
deficit at the position x/c = 4; however we do not observe the vortex that far downstream. The tripping 
of the boundary layer might affect the interaction between the primary and secondary vortices and has been 
reported to decrease the inboard movement of the primary vortex along the span |29| . 

3.5 Temporal evolution 

The Navier-Stokes equations are integrated in time using a second-order accurate stiffly-stable implicit-explicit 
scheme |58| . The advection term is explicitly integrated, whereas the viscous term is implicitly integrated. 
This therefore relaxes the sometimes stringent, diffusion stability condition At oc Ax'^jiy. The CFL restriction 
At cx Ax/u, where u is the advection velocity within each cell, remains however. Because the mesh we consider 
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Figure 5: Contours of iso-helicity showing the interaction between the primary (in light grey) and a secondary 
vortex of opposite rotational sense (in blue, or dark grey). The light grey/blue plane denotes the location of 
the streamwise cut xjc = 0.125 at which results are compared with experiment in Figs. To]|IT 


is coarse, we assume numerical error is dominated by error from the spatial discretization. The time step used 
for the computation normalized by the convective length scale of the chord is At/tc = 1.6 x 10“® which 
translates into a maximum normalized sampling frequency of 3 • 10®. For spectral/ft.p element methods using a 
second order implicit-explicit time integration scheme, the analog to the CFL definition imposes a restriction 
on the maximum timestep of the form|40| At < Aa;/(maxQegQ{|l/®|P^}) where we assume max|F®| ^ Uoo- 
For the present mesh, where the smallest mesh element is 10“^c and using 5*® order polynomials, the timestep 
restriction is of the order of 10“®tc- In practice, since the velocity in some regions can be significantly larger 
than Vaai we use a timestep one order of magnitude smaller. 


4 Results and Discussion 

In this section, the performance of the SVV-iLES method is evaluated by simulating the development of the 
wingtip vortex in the near field and comparing the results against the experimental study by Chow et al. [1] 
as well as the previous LES by Uzun et al. 

We define the near field to be the region above and below the wing and up to one chord length downstream 
of the trailing edge. In this context the mid-field is the region from c to 10c and the far-field is at a distance 
of more than 10c from the trailing edge. To put the challenge of accurately computing the vortex into 
perspective, we outline the typical tracking distances of interest for different applications. The automotive 
industry is interested in tracking the vortex for roughly 20c. For wind turbines and helicopters blades, with 
an aspect ratio of roughly 10, the study of the interaction between the rotor blades and the preceding blades 
that leads to noise and structural vibration, requires tracking the wingtip vortex over more than 60c for one 
revolution. In this study an effort is made to characterize the performance of the modeling method both in 
the three-dimensional turbulent boundary layer and in the rollup wake. It is supposed that accurate modeling 
of the vortex in the near field downstream of the trailing edge pre-supposes an accurate modeling of the 
three-dimensional boundary layer roll-up on the wing surface. 

The vortex can loosely be defined as the region in which the fluid has high helicity, low relative pressure 
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and an axial velocity surplus. Different methods have been developed for defining the center of a vortex, but 
these can prove challenging to apply in the near wake where the vortex is forming as a consequence of the 
shear layer roll-up. These methods include helicity peak correction, vorticity peak, Q'Criterion, zero in-plane 
velocity and axial velocity peak. It has been shown that the helicity peak correction method is best suited 
for estimating the vortex core radius, axial velocity peak and swirl velocity peak. Giuni & Benard [23] assert 
that the centering method chosen should depend on the key aspect of interest but the robustness of these 
methods is not sufficient for identifying the vortex as it develops from the shear layer roll-up and interacts 
with secondary structures. For this reason the method of the manual location of the vortex core by identifying 
local pressure minima has been favored. This method has also been used by Chow et al. so this source of 
error should be taken into account when appreciating the discrepancy between reported core locations (Fig. 
in this region. 

The flow over the wingtip develops into a highly skewed three-dimensional boundary layer that rolls up 
and detaches into a rapidly rotating vortex, at a distance of around 0.5c. This forms an increasingly low 
pressure region in the vortex core that gradually accelerates the fluid entering the core into a jet characterized 
by a notable normalized axial velocity surplus. This strong vortex is thought to be laminar and persistent, 
extending many chord lengths downstream of the trailing edge. The challenge from a modeling perspective 
comes from the three-dimensional boundary layer, the detachment and the strong curvature induced both by 
the geometry and by the many interacting vortical structures. Two key regions of the flow are used to assess 
the performance of the SVV-iLES method against the experimental data of Chow et al. [1]: the wing surface 
and the vortex core. These two regions are of particular interest because they radically differ in nature. Most 
classical turbulence models can capture turbulent boundary layers well, but few of these can accurately model 
the vortex core where curvature of the streamlines is high |59| . 

In the second region of interest, the vortex core, the key feature to reproduce is the low pressure region 
driving the jetting phenomena or the normalized axial velocity surplus. Inside this region, the high value of 
Wmax/Goo at the trailing edge, and the low decay within the first chord length downstream of the trailing 
edge at xjc = 0.867 is a challenging feature to capture. We will therefore compare our results, as well as 
those from previous computations, against the reported experimental values of Umax = l.TTCfyo and 1.59C/oo 
respectively from Chow et. al. [T]. Uzun et al. [^ underpredict the peak normalized and time averaged axial 
velocity by more than 20%, albeit at a lower chord Reynolds number Rcc = 5 • 10®. The k — oj based RANS 
computation by Churchfield et al. |15| report a 150% error, with respect to the experimental data by Chow 
et. al. |T], in the Cp distribution within the vortex at xjc = 0.867. Even the k — to SST-RC RANS [IS], that 
predicts the at x/c = 0.867 with less than 10% error, has nearly 20% error for the estimation of the axial 
velocity surplus and 30% error for the static pressure in the core at cc/c = 0.867. 

The results from SVV-iLES and evaluation with respect to experiment for the two regions of interest is 
presented in the next section. All results presented here have been time-averaged over three chord convective 
length or Itc. We assume the flow is fully developed when the Cp distribution on the wing-surface at the 
span-wise location z/c = 0.899 converges to a smooth curve. 

4.1 Wing Surface 

4.1.1 Cp distribution 

Chow et al. [T| reported the pressure distribution at two spanwise locations. The first cut, at zjh — 0.833, 
is situated inboard of the vortex core where its influence is mild, whereas the second is located at the vortex 
core in zjh = 0.899. At this position, the presence of the vortex leads to a distinctive pronounced suction 
region. The extent of this region upstream depends on the shape of the roll-up layer. The presence of the 
vortex reduces the pressure in this region which translates into an increase in lift. 

Despite a relatively good agreement with experiment for the Cp distribution at the spanwise location 
zjc = 0.833 (Fig. §i). at the spanwise location of the developing primary vortex (Fig. |§)), the SVV-iLES 
computed results under-predict the vortex suction from xjc = 0.5 to xjc = 0.9 and significantly over-predict 
the suction for the last 0.1c. Coarse tangential resolution, in the first 0.5c, may have led to a strong SVV 
dissipation that which in turn significantly damped the early growth of the vortex over the wing surface. The 
sudden change in trend at xjc = 0.9 might be due interaction between the primary and a secondary the 
secondary vortex. With this exception however, the main features of the flow are well captured. 

It should be noted that the instantaneous field is noisy because of the unsteady nature of the boundary layer 
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Figure 6: Comparison with experiment [T] of time-averaged Cp distribution over 3tc time units as a function 
of streamwise position, with the leading edge in a; = 0, at spanwise location zjh = 0.833 (a) and zjh = 0.899 in 
(b) where we also show the results from the uniform grid-refinement study (convergence in P). The number of 
mesh degrees of freedom for 4*'', O*** and 7*^*' order accurate in space are 5.7M, 16.7M and 25.3M respectively. 
The 50% increase in number of degrees of freedom when using 7*^ instead of 6*^ does not significantly affect 
the pressure distribution on the wing at the spanwise location zjc = 0.899. 


at this high Reynolds number and is obviously reduced as we average over longer time intervals, which explains 
the residual noise in the Cp distribution (Fig. ^-b). The use of the spherigon mesh smoothing technique in the 
representation of the geometry may also lead to some higher frequency oscillations and thus less smoothness 
in the distribution. Adding further diffusion from the SVV may help in removing some oscillations. However, 
additional diffusion might lead to an artificial reduction of the Reynolds number. 

4.1.2 Resolution study 

Although this test case is computationally expensive to simulate, we have performed a limited p-refinement 
study of the flow physics, using the Cp distribution as a benchmark for observing convergence and providing 
a form of self-validation. In these tests, the polynomial order was varied, comparing the P = 5 results that we 
present here to results at P = 3 and P = 6. The resulting Cp distributions, presented in figure [^, show very 
little difference between the P = 5 and P = 6 cases, despite a 50% increase in the total number of degrees 
of freedom. However, there is a significant difference between that of P = 3 and P = 5. Whilst there is still 
variation between the experimental results and both P = 5 and P = 6, we can at least conclude that in terms 
of polynomial order the simulation is well-resolved. 

We note that this type of refinement is good in that it is hierarchical and so all the degrees of freedom at 
lower polynomial orders are contained within the higher order simulations. It does not, however, guarantee 
that there cannot be regions which are still not captured and so the basic features of the flow will remain 
reasonably similar. Further work is therefore required to generate a larger sequence of meshes to study the 
effect of refinement in terms of element size. However, given the significant undertaking of this type of study 
and the convergence we have obtained in p, we do not consider this here and use the P = 5 results in the 
coming section. 

4.2 Propagation of the vortex core in the near-wake 

In the vortex core we track the normalized time-averaged axial velocity, the static pressure and at the location 
of the vortex both in the spanwise direction (z) and normal to the suction side direction (positive y), as 
shown in figures and Therefore both time-averaged pressure coefficient and axial velocity have not been 
corrected for the error stemming from possible meandering. In this section we also presents an overview of 
the development of the vortex in the near wake. A comparison shown for both streamlines and time-averaged 
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Figure 7: Comparison with experiment ^ and previous LES by Uzun et al. [5] of the progression of (a) axial 
velocity and (b) Cp distribution in the vortex core. The origin, x/c = 0, is taken to be the position of the 
wingtip trailing edge. 


normalized axial velocity at two crossflow locations: over the wing towards the trailing edge at x/c = —0.114 
and in the near wake at the x/c = 0.125. These results are shown in figures and 

4.2.1 Normalized axial velocity in the vortex core. 

The progression of the axial velocity of the vortex can also be affected by the presence/absence of a boundary 
layer trip as reported by McAlister |55]: for Rcc = 1.5-10® a NACA 0015 profile at angle of attack a = 12° the 
streamwise component of the velocity in the vortex core has a small jetting behavior (< 5% velocity excess) 
without the trip and a significant 20% deficit {u/Uoo ~ 0.8) when a boundary layer trip is added to the 
leading edge. In Fig. ?? we report the streamwise progression of the normalized axial velocity. Our modeling 
adequately resolves the strong jetting behavior measured experimentally. We do however over-predict the 
peak axial velocity in the vortex core by 6% with respect to the experimental value. 

4.2.2 Static pressure within the vortex core. 

In Fig. ?? we report the streamwise progression of the vortex core static pressure. Despite a relative error of 
less than 10% over the wing surface the error grows linearly downstream of the trailing edge reaching 30%. 
Accurately capturing the low pressure region within the vortex core and sustaining this low pressure even just 
one chord distance downstream of the trailing edge is particularly challenging. The linear increase in pressure 
is most likely the result of either too coarse a mesh within the vortex core and/or too strong a contribution 
from the SVV. When the grid is too coarse to adequately capture a gradient the SVV Alter damps out part 
of its kinetic energy. Neglecting compressibility effects, this resulting decrease in jetting velocity increases the 
pressure (or decreases the suction of the vortex). 

4.2.3 Vertical position of the vortex core. 

The vertical location of the primary vortex core is computed to be 20% lower than both the experiment and 
previous LES (Fig.[^). The change in attitude as the vortex leaves the proximity of the wing surface follows 
a similar trend, with the core remaining at the same distance above the wing surface from streamwise location 
xjc = —0.4 from the trailing edge, to the trailing edge xjc = 0 and then showing a pronounced upward trend 
from the leading edge to the experimentally reported xjc = 0.4 location downstream of the leading edge. 
Fig. lit shows the computed cross-section of the vertical vortex profile above the wing surface, with the same 
position measured experimentally (fig. [^) and in the previous LES result of Uzun et al. (flg[^). Despite 
a discrepancy regarding the vertical position of vortex core, we seem to be qualitatively agreeing with the 
experimental results. In particular, the shape of our computed vortex is closer to the isotropic round shape of 
the experimentally obtained vortex, particularly when compared to the previous LES result. The topology of 
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the flow in the region 0.03 < yjc < 0.06, with the presence of a secondary vortex, is also qualitatively closer to 
the experimental results than the previous LES. We should note however that the location of the SVV-iLES 
computed structures is different with respect to those measured experimentally. Whilst it is difficult to explain 
this discrepancy without a further series of detailed simulations, one possible explanation is in the tangential 
grid spacing. We note that for computational reasons, this is clearly not as fine as the wall-normal direction, 
and this may therefore play a key role in influencing the detachment location and therefore vertical position 
of the vortex core. It is also interesting to remark how accurate the LES results from Uzun et al. are at 
predicting the vertical position of the vortex core above the wing despite significantly different flow topologies 
(Fig.|^ and Fig. |^). 




Figure 8: Comparison with experiment [T] and previous LES by Uzun et al. of the vertical position of (a) 
the vortex above the wing and (b) spanwise location. The origin, Zcorejc = 0 and ycoreic = 0 is taken to be 
the position of the wingtip trailing edge. The position of the wing is identified with the rectangle in b). 



Figure 9: SVV-iLES computed streamlines and normalized time-averaged axial velocity at the crossflow plane 
xjc = —0.115 downstream of the trailing edge in b) compared against experimental results from Chow et 
aZ.[T], in (a) and previous LES results by by Uzun et al. 0, in (c). 


4.2.4 Spanwise position of the vortex core. 

The spanwise position of the computed vortex core is compared against both the experimental data of Chow 
et aZ.[T] and previous LES results by Uzun et al. Three key features of flow can be assessed by this figure: 
the location of the origin of the primary vortex, the location at the trailing edge and the evolution of the 
vortex as it leaves the vicinity of the surface of the wing. There is good agreement between experimental 
data and both LES regarding the position of the vortex at the trailing edge. There is, however, a significant 
discrepancy regarding the origin of the vortex. Indeed both Chow et al. [T] (circles in Fig. |^) and Uzun 
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Figure 10: SVV-iLES computed streamlines and normalized time-averaged axial velocity at the crossflow plane 
xjc = 0.125 downstream of the trailing edge in b) compared against experimental results from Chow et al. [T], 
in (a) and previous LES results by by Uzun et al. 0, in (c). Fig. i acts as a companion figure to locate the 
x/c = -1-0.125 with respect to the primary and secondary vortices with a three dimensional perspective of the 
flow in this region. Fig. [D complements Fig. in aiding the identification of the secondary vortex with 
respect to streamlines and also the streamwise component of the vorticity. 


et al. |S] (dot-dashed line in Fig. |^) report the origin of the primary vortex in the region of the wing cap 
whereas the present results show the origin to be on the suction side of the wing around mid-chord. McAlister 
& Takahashi [29] as well as Thompson |M] reported the origin of the vortex on the suction side of a NACA 
0015 profile for a similar case. McAlister & Takahashi [29] also report tripping the boundary layer decreases 
the inboard movement of the primary vortex along the span without altering its distance above the wing. This 
may offer possible insight into the difference in spanwise trajectory between Uzun’s LES computed vortex 
which has a more pronounced inboard movement than Chow’s experiment. 

The most interesting feature of the flow that can be analysed with Fig. is the evolution of the vortex 
as it leaves the vicinity of the surface of the wing. The experimental results of Chow et al. show two distinct 
kinks at xjc = 0.125 and then xjc = 0.452. This is particularly evident when comparing against the previous 
LES of Uzun et al., where the progression of the vortex core moves steadily inboard. We believe these kinks are 
the result of the interaction between the primary vortex and a secondary vortex orbiting around it. Evidence 
of this secondary vortex can be seen in both the comparisons of the streamwise-normal streamlines shown in 
Fig. 10 r with the notable presence of a flat spot in the region 0.72 < zjc < 0.75 and 0.01 < y/c < 0.07. Our 
numerical results also show a similar flat spot, albeit in slightly different position, in region 0.72 < zjc < 0.75 
and 0.06 < y/c < 0.1. By plotting the streamwise component of the vorticity vector in this region we can 
correlate these flat spots to a weaker, counter rotating, secondary vortex (Fig. 111). In Fig. [T^ we enlarge 
this region of interest to evidence the weaker, counter rotating vortex in light grey. In these two figures the 
primary vortex appears in dark grey. This secondary, weaker, counter rotating vortex can also be seen when 
visualising the iso-helicity surfaces. In Fig. the primary vortex appears in grey and the secondary vortex in 
dark grey (or blue if visualised in color). The light grey (light blue) surface represents the spanwise location 
x/c = -1-0.125 at which the comparison is made between experiment, previous LES and present results for 
Fig.[TgTT| It is also interesting to note that the computed secondary vortex seems to be out of phase, in 
the streamwise direction, with respect to the experimentally computed vortex. Indeed it appears in the II 
quadrant whereas in the experimental results it appears in the III quadrant (when viewed as in Fig. [Toplt . 


5 Conclusion 

The SVV-iLES workflow, developed for computing unsteady vortex-dominated flows, has been assessed by 
comparing numerical results with experimental data by Chow et aZ.[T] for a NACA 0012 wingtip vortex test 
case, at a higher Reynolds number than any LES study performed to date. Overall, the results show the 
potential of this method to resolve the large scale features of the flow, without the use of explicit turbulence 
and sub-grid scale models. The use of an implicit LES presents a notable advantage over these methods, given 
that only two parameters are needed to control regularization and stability of the numerics. 
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Figure 11: Streamwise component of the vorticity vector at location xjc = +0.125. Both primary vortex, 
in dark grey, and secondary vortex in light grey in the II quadrant are visible with an overview of the flow in 
a) and the detail of the location of the secondary vortex in b). These figure acts as a companion to Fig. 10 d. 


Our results show better correlation with experimental results than previous numerical results, both in 
terms of the static pressure distribution, prediction of the jetting velocity, vortex spanwise location and the 
ability to resolve the secondary vortex interaction with the main wingtip vortex. In particular, we note that 
the results presented here show a good agreement in the roll-up region where both the static pressure and 
velocity magnitude agree within 10% of experiment. We also observe that although the mesh used in this 
study is coarser than ones reported in other studies, as highlighted in Tab. and the Reynolds number is 
larger than previous investigations, the results of this study demonstrate that the SVV-iLES method can still 
accurately capture the essential features of the flow. However we do note that unsteady simulation obviously 
requires significantly more compute resource as compared to steady RANS simulations. We additionally 
obtain a qualitatively good agreement of the modeling of the vortex roll-up, and predict the secondary vortex 
and its interaction with the primary vortex over the wingtip. We believe this can be attributed both to the 
lower diffusion and dispersion properties of the spectral//ip element method and to the use of an isoparametric 
refinement technique which provides adequate wall-normal resolution in the sub-viscous layer of the boundary 
region. 

There are however some clear differences between the results presented here and experimental data. It is 
clear that the suction peak on the wing appears further downstream {xjc = 0.95 instead of the experimental 
value of xjc = 0.85). There is also a visible difference in the location of the secondary vortex at location 
xjc = 0.867 downstream of the trailing edge. These points seem to indicate that despite successfully 
modeling the secondary vortex, the interaction between the primary and secondary vortex is not yet accurate 
enough to reproduce experimental results. The p-reflnement study presented here shows that, whilst our 
results are well-resolved in terms of the polynomial space, a small increase in polynomial order does not 
generally lead to a better correlation with the experimental data. This is likely due an under-resolution of 
the wall-tangential directions across of the surface of the wing and in particular in the first 0.5c, where the 
primary vortex originates. Therefore, whilst a large increase in polynomial order would hopefully yield a 
better correlation with the experimental results, a more efficient approach to achieve convergence is likely 
to be a combination of both mesh and p-refinement. It is therefore clear that, together with improving the 
smoothness of the surface mesh and a further increase in Reynolds number to match the experiment, future 
studies should include additional local refinement in terms of element size in order to hopefully attain a closer 
agreement with the experimental results. 

In summary, the SVV-iLES method has been shown to be a compelling alternative for computing complex 
unsteady vortex dominated flows, such as the wingtip vortex, motivating its use for complex industrially 
relevant cases where high-fidelity computational fluid dynamics can become an enabling technology. 
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